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ABSTRACT 

We study the dependence of the delayed neutrino-heating mechanism for core-collapse supernovae on the 
equation of state. Using a simplified treatment of the neutrino physics with a parameterized neutrino luminosity, 
we explore the relationship between explosion time, mass accretion rate, and neutrino luminosity for a 15 Mq 
progenitor in ID and 2D. We test three different equations of state commonly used in core-collapse simulations: 
the models of Lattimer & Swesty (1991) with incompressibility of 180 MeV and 220 MeV, and the model of 
Shen et al. (1998), in order of increasing stiffness. We find that for a given neutrino luminosity the time after 
bounce until explosion increases with the stiffness of the equation of state: the Lattimer & Swesty EOS explode 
more easily than that of Shen et al. We find this holds in both ID and 2D, while for all models explosions 
are obtained more easily in 2D than in ID. We also discuss the relevance of approximate instability criteria to 
realistic simulations. 

Keywords: supernovae: general - hydrodynamics - neutrinos - stars: interiors 



1. INTRODUCTION 

The exact process that halts the collapse of the core of a mas- 
sive star at the end of its life and drives a successful supernova 
explosion, releasing a multitude of neutrinos and ejecta with 
around 10^^ erg of kinetic energy is not yet fully-understood. 
The delayed neutrino-heating mechanism of core-collapse su- 
pernovae (CCSNe, Colgate & White 1966; Bethe & Wilson 
1985), the leading candidate for the explosion mechanism, 
fails to consistently produce energetic explosions for the wide 
range of possible progenitor masses in simulations that include 
high-fidelity treatments of neutrino transport. Neutrinos are 
an attractive explanation for the explosion mechanism as the 
bulk of the gravitational binding energy of the progenitor core, 
some 10^^ erg, will be radiated away as neutrinos as the proto- 
neutron star (PNS) forms and cools; only about 10^^ erg of this 
energy needs to be transferred to the collapsing stellar material 
in order to explain typical CCSNe energies. Only a fraction, 
however, of the core's original gravitational binding energy, 
perhaps a few times 10^^ erg, will be released as neutrino radia- 
tion in the first second following the collapse and bounce of the 
core, approximately the timescale on which a successful explo- 
sion must occur. Neutrino radiation-hydrodynamic simulations 
of core collapse show that this is challenging to achieve, partic- 
ularly in ID simulations. It is now clear that multi-dimensional 
effects, such as proto-neutron star convection (Epstein 1979; 
Burrows & Fryxell 1993; Dessart et al. 2006; Burrows et al. 
2007a), neutrino-driven convection (Janka & Mueller 1996), 
the standing accretion shock instability (SASI, Blondin et al. 
2003), and turbulence (Murphy & Meakin 201 1) play a critical 
role in the success of the neutrino mechanism in driving explo- 
sions. Together, these multidimensional effects can push some 
progenitors over the critical threshold, resulting in somewhat 
marginal explosions in 2D in certain simulations (Marek & 
Janka 2009; Bruenn et al. 2009; Suwa et al. 2010; Mueller et al. 
2012). This result is dependent on details of the numerical 
scheme and treatment of neutrino transport; the 2D simulations 
of Burrows et al. (2006, 2007b) do not find neutrino-driven 
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explosions for any progenitors.^ Based on a series of general 
relativistic ID simulations including Boltzmann transport for 
neutrinos, Lentz et al. (2012) suggest that the discrepancies in 
the results from the various groups simulating CCSNe may be 
due in large part to differences in how three major aspects of 
the physics are treated: general relativistic versus Newtonian 
dynamics, the neutrino interactions and opacities considered, 
particularly inclusion of inelastic neutrino scattering, and ob- 
server frame corrections in the neutrino transport. 

All of this indicates, unsurprisingly, that simulations of the 
neutrino mechanism are sensitive to microphysics, and details 
of how the microphysics is treated. It is still uncertain, how- 
ever, if the variation in the handling of microphysics is large 
enough to explain the dearth of energetic explosions; that is, if 
a more accurate handling of the details of the neutrino physics 
will result in successful explosions having around 10^ erg of 
kinetic energy for a wide range of progenitor masses. Nord- 
haus et al. (2010) have suggested a larger impact on the success 
of the neutrino mechanism results from fully 3D simulations. 
Using a simplified neutrino heating and cooling scheme rather 
than expensive neutrino transport, Nordhaus et al. (2010) find 
that the driving neutrino luminosity necessary to obtain an 
explosion is about 20% lower in 3D than in 2D. This differ- 
ence may be larger than what we could expect to gain from 
higher-fidelity treatments of neutrino transport or inclusion of 
general relativity in simulations. Hanke et al. (2011) sought to 
reproduce this result. Employing a neutrino heating/cooling 
approach similar to Nordhaus et al. (2010), Hanke et al. (201 1) 
find the critical neutrino luminosity in 3D is not significantly 
lower than that of 2D. This raises the questions of just how 
beneficial 3D simulations are to the success of the neutrino 
mechanism and what is the cause for the differences in the 
results of Nordhaus et al. (2010) and Hanke et al. (2011)? 
As 3D simulations become more feasible and available, we 
may be able to answer the question of the importance of 3D. 

^ Burrows et al. (2006, 2007b) find explosions via the "acoustic" mecha- 
nism first revealed by their simulations. The work of Weinberg & Quataert 
(2008), however, shows that the ^ = 1 g-mode vibration of the neutron star 
which powers the acoustic mechanism will in reality saturate at energies well 
below that required to drive an explosion. 
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Table 1 

EOS Parameters and experimental limits. 



EOS 


(MeV) 


(MeV) 


(MeV) 


(MeV) 


/\/fmax e 

(Mo) 


(km) 


LSI 80 


180 


-451 


28.6 


74 


1.84 


12.2 


LS220 


220 


-411 


28.6 


74 


2.06 


12.7 


LS375 


375 


-162 


28.6 


74 


2.72 


13.5 


STOS 


281 


-285 


36.9 


111 


2.22 


14.6 


HS (TMA) 


318 


-572 


30.7 


90 


2.02 


13.5 


Limits 


240 ± 10 


-355 ± 95 


~ 32 


~ 75 


> 1.97 ±0.04 


11 - 12 



^ Nuclear incompressibility. 
^ Skewness coefficient. 
^ Symmetry energy coefficient. 
^ Symmetry energy slope. 

^ Maximum gravitational mass of a cold neutron star 
^ Radius of a 1.4 Mq neutron star 




Figure 1. Mass-Radius relations for cold neutron stars as computed by solving 
the Tolman-Oppenheimer-Volkov equations for the various EOS Hsted. 

Recently, Takiwaki et al. (2012) report an explosion for an 
11.2 Mq progenitor in a 3D simulation performed using the 
isotropic diffusion source approximation for neutrinos (IDSA, 
Liebendorfer et al. 2009). Compared to 2D, Takiwaki et al. 
found that the convection below the gain region was more 
vigorous in 3D enhancing the neutrino luminosity beyond the 
2D case. The somewhat modest resolution they were able to 
afford, however, means that the issue of the importance of 3D 
simulations cannot be settled on the basis of their results alone. 

Besides the neutrino physics, an important microphysical de- 
pendence of the CCSN mechanism is the equation of state. The 
equation of state for matter at densities and temperatures rele- 
vant to CCSNe is still an active area of research (e.g., Hempel 
& Schaffner-Bielich 2010; Shen et al. 2010, 2011a,b,c; Furu- 
sawa et al. 201 1). Our understanding of nuclear interactions at 
these densities and temperatures is incomplete and a number of 
models exist. Testing the various EOS models with experiment 
and observation is challenging, but a few constraints do exist. 
Hempel et al. (2012) review some of the constraints on the 
EOS parameters from laboratory experiments. In Table 1 we 
summarize the constraints on the nuclear EOS parameters and 
the respective numbers for five EOS models used in CCSN 
simulations: Lattimer & Swesty with incompressibility param- 
eters of 180 MeV (LSI 80), 220 MeV (LS220), and 375 MeV 
(LS375), Shen et al. (STOS), and the recent EOS of Hempel 



& Schaffner-Bielich (2010) using the TMA model for nuclear 
interactions (HS (TMA)). Most values for the limits are taken 
from experiment (for references to the experimental data see 
Hempel et al. 2012), except for the skewness coefficient, K\ 
which is only a theoretical estimate. Based on the measure- 
ments listed in the table, only LS220 falls within or near the 
limits for each of the EOS parameters. 

Observations of neutron stars and pulsars can also place 
constraints on the EOS for hot dense matter. Using a parame- 
terized EOS, Steiner et al. (2010) analyze a small sample of 
neutron stars and determine a most likely mass-radius rela- 
tionship for cold neutron stars. Their results favor a radius of 
11 — 12 km for a 1 .4 Mq neutron star, and a maximum neutron 
star mass in the range 1.9 — 2.2 Mq. The maximum mass 
of a cold neutron star is an important constraint for the EOS 
models. The observations of PSR Jl 6 14-2230 by Demorest 
et al. (2010) place a very tight limit on the mass of the neu- 
tron star in this system of 1.97 ± 0.04 Mq. Given an EOS 
model, solutions to the Tolman-Oppenheimer-Volkov equation 
yield the gravitational mass-radius relationship for cold neu- 
tron stars. Figure 1 shows the mass-radius relationships for 
the EOS of Lattimer & Swesty (1991) and Shen et al. (1998), 
along with the mass limits for PSR J1614-2230. Based on 
the measurement of Demorest et al. (2010), LSI 80 is ruled 
out as it does not produce a large enough maximally-massed 
neutron star. The Shen et al. EOS produces a maximum neu- 
tron star mass well above the limit placed by PSR Jl 6 14-2230 
(2.22 Mq), but STOS predicts a radius for a 1.4 NS of 
around 15 km, well outside the limits estimated by Steiner 
et al. (2010). The theoretical predictions for the radius of a 
1.4 Mq NS of Hebeler et al. (2010) are somewhat broader: 
9.7 - 13.9 km, but STOS still Hes outside of this range. Only 
LS220 satisfies the constraint on the maximum NS mass from 
PSR J1614-2230 and produces a 1.4 M© NS with a radius 
(~ 12.7 km) close to the upper estimate of Steiner et al. (2010) 
and well within the limits of Hebeler et al. (2010). Figure 1 
also shows the mass-radius relationship for LS375. This EOS 
easily satisfies the requirement on the maximum neutron star 
mass, but falls above the radius limits of Steiner et al. (2010) 
and has an incompressibility (375 MeV) well above that found 
in experiment. 

In this work, we study the dependence of the neutrino mech- 
anism on the equation of state employed. Using ID and 2D 
hydrodynamic simulations with simplified neutrino transport, 
we determine the explosion times for a 15 Mq progenitor as 
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a function of neutrino luminosity for three EOS: Lattimer & 
Swesty (1991) with K = 180 MeV and K = 220 MeV, and 
Shen, Toki, Oyamatsu, & Sumiyoshi (1998). In ID only, we 
also run simulations using LS375. Our treatment of the neu- 
trino physics is similar to those used by Murphy & Burrows 
(2008), Nordhaus et al. (2010), and Hanke et al. (2011): we 
assume local neutrino heating and cooling, based on the rates 
derived by Janka (2001), and a constant neutrino luminosity. 
We find that the explosion times for a given neutrino lumi- 
nosity are significantly dependent on the EOS used, with the 
general trend of easier explosions for softer equations of state. 
Thus, the Lattimer & Swesty models result in easier explo- 
sions than the Shen et al. EOS, with LSI 80 leading to the 
earliest explosions at a given neutrino luminosity. Early ID 
simulations using a simple, parameterized EOS showed that 
stronger shocks result for softer EOS (Baron et al. 1985a,b). 
Our results show that this trend is recovered for more realistic 
nuclear EOS and in 2D simulations. 

The EOS has been shown to influence the results of CCSN 
simulations in a number of contexts. Thompson et al. (2003) 
explored the sensitive of their ID neutrino radiation hydrody- 
namic simulations to the value of K for the Lattimer & Swesty 
EOS. They found that the overall evolution during the first 
200 ms post-bounce is not very dependent on K: the resulting 
temperatures were slightly higher and the emergent neutrino lu- 
minosities at most 9% larger for LSI 80 as compared to LS375. 
Sumiyoshi et al. (2005) compared LSI 80 and STOS in ID 
general relativistic simulations. While neither model explodes 
in their simulations, they find differences in composition re- 
sulting from the different symmetry energies (see J in Table 1) 
and in temperature evolution with LS 1 80 giving slightly larger 
temperatures and, therefore, slightly higher neutrino luminosi- 
ties. The higher neutrino luminosity and heating rate found 
for LSI 80 is somewhat counteracted by the more compact 
proto-neutron star that results from the lower incompressibility. 
Sumiyoshi et al. (2005) find that the differences between the 
two models grow with time, particularly after 300 ms post- 
bounce. Marek et al. (2009) explored the dependence on the 
gravitational wave signal produced by 2D simulations on the 
EOS, comparing LSI 80 and the stiffer EOS of Hillebrandt & 
Wolff (1985). They find that the matter-generated gravitational 
wave signal is not very dependent on the EOS but that the 
differences in the anisotropic neutrino emission gravitational 
wave signal, which can dominate the matter gravitational wave 
signal, are significantly different between LS 1 80 and Hille- 
brandt & Wolff EOS. In failed explosions, the formation of 
black holes has been found to be highly dependent on the EOS 
(Sumiyoshi et al. 2007; Fischer et al. 2009; O'Connor & Ott 
201 1 ; Hempel et al. 2012). Generally, stiffer EOS result in a 
greater time to black hole formation following bounce. 

This paper is organized as follows. In Section 2 we describe 
the details of our numerical approach. In Section 3 we describe 
the results of our simulations and the dependence of the explo- 
sions on the EOS. In Section 4 we conclude and discuss how 
our work relates to recent work on neutrino- driven supernovae. 



2. COMPUTATIONAL METHODOLOGY 

We conduct our study using the FLASH^ simulation frame- 
work (Dubey et al. 2009) to solve the Eulerian equations of 
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Figure 2. Comparison of the fitting formula used in this work to compute the 
neutrino optical depth and the integrated optical depth from a one-dimensional 
simulation at 200 ms post-bounce. The black curve is the exponential of the 
negative neutrino optical depth as computed using the effective of opacity of 
Janka (2001). The red curve is the same but using the fitting formula. Plotted 
also is an exponential cutoff based on the density in units of 10^^ g cm~^, 
pii . The bottom panel shows the error between the approximate approaches 
and that based on the computed optical depth. The fitting formula used in 
this work matches the integrated optical depth well, better than the simple 
density-based approach. 



hydrodynamics, 



dp 

dt 



■v-(H = o, 



dpv 

~dt 
dpE 



(1) 

V • {(mv) + VP = -/jV$, (2) 
■V-[{pE + P)v] = pv-'^^ + p{'H-C), (3) 



' Available for download at www.flash.uchicago.edu. 



where p is the mass density, v the velocity vector, P the pres- 
sure, <l> the gravitational potential, E the total specific en- 
ergy, and — C) is the difference between neutrino heating 
and cooling. We use the directionally-unsplit hydrodynamics 
solver in FLASH to solve equations (1) - (3) (Lee & Couch, 
in prep, see also Lee & Deane 2009). The unsplit solver in 
FLASH has many options for spatial reconstruction scheme, 
Riemann solver, slope limiter, etc. For this work we use third- 
order piecewise-parabolic spatial reconstruction (PPM, Colella 
& Woodward 1984), a "hybrid" slope hmiter, and the HLLC 
Riemann solver. Inside of shocks, we use the HLLE Riemann 
solver in order to avoid the odd-even decoupling instability 
in 2D (Quirk 1994). This approach allows us to avoid the 
smearing of contact discontinuities that the HLLE solver is 
prone to while benefiting from the HLLE solver's immunity 
to odd-even decoupling at shocks. The "hybrid" slope limiter 
(Balsara 2004) applies the monotonized central (mc) limiter to 
linear wave families (such as density) and the robust, slightly 
more diffusive minmod limiter to non-linear, self- steepening 
wave families. The gravitational potential, <l>, is calculated 
using the new multipole Poisson solver available in FLASH4 
(see section 8.10.2.2 of the FLASH4 User's Guide), however, 
we set ^max = yielding monopole gravity. 

Following the approaches of Murphy & Burrows (2008), 
Nordhaus et al. (2010), and Hanke et al. (2011) we use the 
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Figure 3. Mass accretion rate as a function of time relative to bounce for the 
15 Mq progenitor used in this study. The mass accretion rate is measured at a 
radius of 500 km for a non-exploding model. The accretion rates for the three 
EOS used are shown. The Lattimer & Swesty EOS differ from that of Shen et 
al. only at early times and are essentially identical after 120 ms post-bounce. 



neutrino heating and cooling rates derived by Janka (2001), 
n = 1.544 X 10^° ' 



1052 ergs-1 ) V4MeV ) 



100 km 



and 

C = 1.399x10^° 



2MeV 



erg 
g • 



erg 
g • 



(4) 



, (5) 



where Lj^^ is the electron neutrino luminosity (note it is as- 
sumed that Lj)^ = LyJ, Tjy^ is the electron neutrino temper- 
ature, r is the spherical radius, (Yp + Fn) is the sum of the 
neutron and proton number fractions, r^^ is the electron neu- 
trino optical depth, and T is the matter temperature. In all 
of our calculations, Tj^^ is set to 4 MeV. The heating term, 
equation (4), is applied only after core bounce has occurred 
in our simulations. The factor e~^''^ is included to provide an 
exponential cutoff of the heating and cooling at high density 
where the material is optically-thick to neutrinos. As such, we 
do not calculate the optical depth self-consistently but instead 
rely on a piece- wise fit of the optical depth as a function of 
density. In Figure 2 we compare our piece-wise fitting ap- 
proach for Tjy^ with the integrated electron neutrino optical 
depth computed using the effective opacity given by Janka 
(2001), A^eff = 1.5 X 10-^10 (A:T^e/4 MeV)^. As shown, 
the piece- wise fit closely approximates the actual integrated 
optical depth. We also show in Figure 2 an exponential cutoff 
based directly on the density. We see that our piecewise ap- 
proach better approximates the integrated optical depth than 
this simple density based approach. We stress that the factor 
e~^^e in Equations (4) and (5) is used simply as a cutoff for 
the heating and cooling functions at high-density and the exact 
values of the cutoff are less important than consistent applica- 
tion and use in all models considered in a study such as the 
present. 

We calculate the evolution of the electron fraction, Y^, ac- 
cording to the parameterization of Liebendorfer (2005) in 



Table 2 

Time post-bounce and accretion rates at time of explosion for the 15 Mq 
progenitor with 0.5 km resolution. 



LSI 80 



LS220 



STOS 





f b 




^exp 


Mexp 


^exp 


Mexp 


(10^2 gj-g/g^) 


(ms) 


(Mo/s) 


(ms) 


(Mo/s) 


(ms) 


(Mq/s) 




ID 


2.0 














2.1 














2.2 


689 


0.192 


816 


0.169 






2.3 


374 


0.254 


571 


0.211 


943 


0.153 


2.4 


213 


0.289 


402 


0.251 


829 


0.170 


2.5 


193 


0.312 


196 


0.311 


380 


0.262 


2.7 


175 


0.317 


178 


0.320 


216 


0.310 


2.9 










200 


0.314 


2D 


1.3 


813 


0.182 










1.5 


277 


0.279 


368 


0.262 


732 


0.192 


1.7 


203 


0.307 


246 


0.295 


374 


0.263 


1.8 






194 


0.313 






1.9 


178 


0.315 






249 


0.298 


2.0 


179 


0.318 


171 


0.322 


217 


0.313 


2.1 






171 


0.322 


201 


0.317 


2.2 






167 


0.324 


193 


0.322 


2.4 










183 


0.322 


2.6 










168 


0.328 



^ Electron-neutrino luminosity. 

^ Time after bounce of onset of explosion. A symbol indicates that the model 
does not explode during the simulated period of evolution. 
^ Mass accretion rate at onset of explosion. 



which Ye is dependent only on density. This approach is ap- 
propriate only for the pre-bounce collapse phase, however 
we employ the parameterization post-bounce as well in order 
to maintain consistency with Nordhaus et al. (2010). With 
no post-bounce evolution of the electron fraction the cooling 
above the neutrinosphere is grossly underestimated leading to 
significantly earlier and easier explosions, as mentioned by 
Hanke et al. (201 1). We do not include the entropy increases 
due to neutrino thermalization as prescribed by Liebendorfer 
(2005). 

We run our simulations in ID spherical and 2D cylindrical 
geometries using adaptive mesh refinement (AMR) as imple- 
mented in FLASH via PARAMESH (v.4-dev, MacNeice et al. 
2000). For this work we have extended FLASH'S unsplit 
solver to spherical and cylindrical coordinates using an ap- 
proach similar to those described in Skinner & Ostriker (2010) 
and Appendix A of Mignone et al. (2005). Our fiducial set 
of simulations has a minimum grid- spacing of 0.5 km in each 
direction, obtained using 8 AMR levels. The grid in our ID 
spherical simulations covers radii from km to 5000 km. In 
2D cylindrical runs, we simulate a grid that spans 5000 km in 
R and 10,000 km in z. We do not remove the proto-neutron 
star from our grid and follow the evolution from progenitor 
phase through collapse and bounce and into the shock stag- 
nation/revival phase. For boundary conditions at the outer 
edge of the domain, we assume radial power-law profiles that 
approximate the stellar envelope. We find that simple "outflow" 
boundary conditions, which enforce a zero-gradient condition, 
are not appropriate for the core-collapse context because the 
amount of stellar material entering the outer edge of the do- 
main will be overestimated, artificially enhancing the mass 
accretion rate at late times. By experimentation we have found 
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Figure 4. Neutrino luminosity versus mass accretion rate at time of explosion (left panel) and time of explosion (right panel) for the three EOS considered in this 
work. Results from both ID (dashed lines) and 2D (solid lines) are shown. For each EOS explosions are found more easily in 2D than in ID. The Lattimer & 
Swesty EOS result in easier explosions in both ID and 2D than the Shen et al. EOS. 




tph (s) V (s) 



Figure 5. Shock radii as a function of time for STOS, LS220, and LSI 80 in ID (left) and 2D (right). Three different neutrino luminosities are plotted for each 
EOS in each panel, as labeled. 



that "outflow" boundary conditions can, in this way, suppress 
explosions for neutrino luminosities near critical. For all of 
our simulations, we use the 15 Mq progenitor of Woosley & 
Weaver (1995). 

We have incorporated the finite temperature equation of state 
routines of O'Connor & Ott (2010) into the FLASH frame- 
work"^. We use three different EOS models in our simulations, 
the models of Lattimer & Swesty (1991) with incompressibil- 
ity, K, of 180 MeV and 220 MeV and that of Shen et al. (1998). 
Some parameters for these EOS are listed in Table 1 . 

3. RESULTS 

3.1. Dependence on EOS 

We ran a series of core collapse simulations in ID and 2D in 
which we varied the driving neutrino luminosity. For models 
that explode, we measure the post-bounce time of the explo- 
sion and the mass accretion rate at the time of explosion. We 
consider a model to have exploded if the average shock radius 
exceeds 400 km and does not subsequently fall back below 
400 km. We measure the mass accretion rate at a radius of 

^ These routines are available for download at stellarcollapse.org. 



500 km in both ID and 2D simulations, though the measured 
mass accretion rates are identical since the 2D solution remains 
spherically-symmetric outside of the shock. Figure 3 shows 
our measured mass accretion rates as a function of time post- 
bounce for the three EOS we consider. Our mass accretion rate 
is very similar to that of Hanke et al. (201 1) for their 15 Mq 
progenitor model. The differences in mass accretion history 
between STOS and LS can be attributed to models using STOS 
collapsing and bouncing a little faster than models using LS; 
bounce occurs 50 ms earlier for STOS as compared to LS. 
Table 2 summarizes the simulations we ran and the resulting 
explosion times and mass accretion rates at the time of explo- 
sion. In Figure 4 we plot the driving neutrino luminosities as a 
function of explosion time and mass accretion rate at the time 
of explosion. 

Our results show that the Lattimer & Swesty EOS explode 
more easily than that of Shen et al., with LSI 80 resulting in 
the earliest explosions for a given neutrino luminosity (lowest 
curves in Fig. 4). The results then follow a basic trend that the 
stiff er the EOS, the harder it is to drive an explosion (LSI 80 
being the softest EOS and STOS being the stiffest EOS we 
consider). This trend holds in both ID and 2D simulations. For 
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Figure 6. Radial entropy plots comparing ID (red lines) and 2D (black dots) for LS220 (left panels) and STOS (right panels) at two different times post-bounce: 
50 ms (top) and 200 ms (bottom). The neutrino luminosity for all the models shown in this figure is 2.0 x 10^^ erg s"-*^. 



high neutrino luminosities, resulting in earlier explosions, the 
curves in Figure 4 for LSI 80 and LS220 converge, showing 
very little difference in time of explosion at these luminosi- 
ties. Figure 5 shows the (average) shock radii in ID (2D) 
for the three EOS and several neutrino luminosities. The ID 
simulations show a clear delineation of models into the three 
categories of the radial shock instability (Fernandez 2012): os- 
cillatory/stable^ , oscillatory/unstable, non-oscillatory/unstable. 
In 2D, £ > modes of the SASI dominate the purely ra- 
dial (£ = 0) shock oscillations seen in ID leading to a more 
chaotic variation in the average shock radius in time. In Fig- 
ure 6 we compare the radial entropy profiles in ID and 2D 
for LS220 and STOS at 50 ms and 200 ms post-bounce for 
Lu^,52 = ^zye/(10^^ erg s~^) = 2.0. For this luminosity, both 
2D simulations result in rather prompt explosions whereas 
the ID simulations fail to explode entirely with the shock 
stalling around 200 km. The entropy scatter plots for the 2D 

^ In our simulations, since the neutrino luminosity remains constant while 
the mass accretion rate decreases with time, there is really no such thing as a 
truly stable oscillatory model. All models at a given neutrino luminosity will 
eventually explode once the mass accretion rate has dropped sufficiently. 



simulations shown in Figure 6 also demonstrate the entropy 
re-distribution due to convection and turbulence. 

Figures 7-9 show the evolution of the entropy in 2D for 
Liy^ = 1.3 X 10^^ erg s~^. For this luminosity, only LSI 80 
results in an explosion within 1 second post-bounce. Figure 7 
corresponds to 100 ms post-bounce. At this stage convection 
has set in in the post-shock region as well as the PNS cooling 
region (faint blue area in the figure). Some extension along 
the axis is evident, particularly for LS 1 80, but a strong £ = 1 
SASI is not yet apparent. By 300 ms (Figure 8), however, 
the SASI is very obvious and is reaching large amplitude for 
LSI 80, though less so for LS220 and STOS. At 600 ms (Figure 
9) large-amplitude £ = 1 SASI motions are clear for LS 1 80 
and LS220, but the SASI amplitude remains modest for STOS. 
While LS 1 80 and LS220 are showing shock expansion, partic- 
ularly along the axis, STOS shows a fairly stationary shock, in 
an average sense, at around 200 km. 

Our results show a clear trend with the incompressibility, K, 
of the EOS: higher values of K result in later explosions. In 
order to test the importance of K in determining the time to 
explosion, we have also run a series of ID simulations using 
the Lattimer & Swesty EOS with K = 375 MeV. The results 
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are shown in Figure 10 where we plot the neutrino luminosity 
versus mass accretion rate at time of explosion for our ID 
simulations. As can be seen, the trend with K does not survive: 
the simulations with LS375 explode earlier (at higher M) than 
STOS. At low-luminosities, near critical, LS375 explodes later 
than LS220 and LSI 80. At higher luminosities, LS375 roughly 
converges with LS220, and all three versions of the LS EOS 
converge at Ly^^^2 — 2.5. This indicates that the time to 
explosion for a given EOS is not dependent on K alone, but 
must also depend on the higher-order EOS parameters listed 
in Table 1 . 

Why, given an identical neutrino luminosity, do differences 
in the baryonic EOS result in different evolution of the super- 
nova shock? That is, what mechanism couples the behavior of 
the PNS, where the differences between the EOS are important, 
to the motion of the shock? The difference in time-to-explosion 
we find for the various EOS is likely due to EOS-dependent 
behavior of the SASI. Whether the underlying SASI mech- 
anism is solely acoustic (Blondin & Mezzacappa 2006) or 
advective-acoustic (Foglizzo et al. 2007) the magnitude of the 
acoustic energy flux from the PNS is a determining factor in 
its growth (Guilet & Foglizzo 2012). ^ A soft EOS, such as 
LSI 80 and LS220, results in a more compact, higher-density 
PNS as compared to a stiff EOS like STOS. This is true for the 
central core of the PNS as well as the 'deceleration' region at 
the edge of the PNS where advected perturbations are turned 
into sound waves. The EOS dependence we find can likely 
be explained by a difference in the acoustic flux generated at 
the edge of the PNS: the denser PNS's produced by LSI 80 
and LS220 generate a larger outbound acoustic flux than that 
of STOS and, thus, lead to larger amplitude excitations of the 
SASI. Dense PNS are, in effect, better sonic transducers. This 
will hold in both 2D and ID, where the SASI is purely radial. 
Differences in the growth of the SASI also accounts for the 
similarity in time-to-explosion for LSI 80 and LS220 at high 
Ly : these prompt explosions do not depend so much on the 
SASI. 

The higher mass accretion rate of STOS at early times (Fig- 
ure 3) is causal to the differences in time-to-explosion for 
high-L^^ simulations. This higher accretion rate prior to 100 
ms post-bounce results in a greater mass of iron being disso- 
ciated at the shock and, therefore, more energy loss. In the 
high-L^^^ models, the shock begins running away to large radii 
at the moment the mass accretion rate drops precipitously (see 
Figure 5), corresponding to the end of the accretion of the iron 
core. Since this is delayed slightly for STOS, the shock be- 
gins its runaway expansion later as compared with LS 1 80 and 
LS220. This difference will be less important at later times and 
cannot explain the differences in time-to-explosion between 
LSI 80 and LS220 (and LS375 in ID) we find for lower-L^^ 
models. 

3.2. Criterion for Transition to Runaway Shock Expansion 

Recently, Pejcha & Thompson (2012) and Fernandez (2012) 
have studied the conditions that lead to explosions in simple, 
parameterized models of stalled shocks in the CCSN context. 
Using principles of stellar pulsation theory, Fernandez (2012) 
develops an instability criterion that closely resembles the 
critical condition found in previous simulations: a greater ad- 
vection time through the gain region than through the cooling 

^ Guilet & Foglizzo (2012) present a compelling argument that the linear 
growth of the SASI must be dominated by the advective-acoustic cycle rather 
than the purely acoustic cycle. 
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Figure 10. The neutrino luminosity versus mass accretion rate at time of 
explosion for ID simulations using four different EOS, including LS375. 

region. Pejcha & Thompson (2012) find a criterion based on 
the ratio of the sound speed to the escape velocity, the so-called 
'ante-sonic condition.' We find that both criteria can be good 
predictors for the transition to runaway shock expansion in our 
ID simulations. The ante-sonic condition is, however, easier 
to compute (especially in 2D) since it does not involve spatial 
integrals. The ante-sonic condition for a polytropic EOS is 

max {^^^ = O-lSr, (6) 

where ^esc is the escape velocity, and V is the polytropic in- 
dex. In Figure 11 we plot the quantity (c|/Vesc)/(0-lS7c) as 
a function of space and time for the ID simulation using the 
STOS EOS with Ljy^^^2 = 2.3, along with the shock radius as 
a function of time for this model. For the general EOS, we 
substitute the spatially- varying adiabatic index as returned by 
the EOS, 7c, for the polytropic index in the ante-sonic condi- 
tion. In Figure 11, instability of the type described by Pejcha 
& Thompson (2012) occurs for values greater than one. The 
highest values for the ante-sonic relation occur when the shock 
reaches radial maxima, thus we plot the times when the shock 
has just started to collapse back during its radial oscillations. 
As shown, the ante- sonic condition is only satisfied on the last, 
explosive expansion of the shock. This is true for all of our 
ID models except the Ly^^^2 =2.4 model using STOS; in this 
model, we see the shock expand well beyond 400 km and then 
fall back inward one last time before runaway shock expan- 
sion occurs (see Figure 5). During the shock's first expansion 
beyond 400 km in this model, the ante- sonic condition is met 
but the shock nevertheless contracts again. 

In 2D, we find that the ante- sonic condition is generally 
met in post-shock regions that have expanded to large radii. 
This can happen long before the average shock radius exceeds 
400 km when we consider the model to have successfully 'ex- 
ploded.' For radial rays along which the ante-sonic condition 
has been met at some point, the shock radius tends to remain 
above 400 km. 

3.3. Dependence on Resolution 

Hanke et al. (201 1) report that their simulations, which em- 
ploy the same parametric neutrino heating/cooling as this study, 
are significantly resolution-dependent. In their 2D simulations, 
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Figure 11. The ante-sonic condition as a function of radius at several times post-bounce for the STOS model with L^y^ = 2.3 x 10^^ erg s ^ (left panel). Shown 
for comparison is the shock radius as a function of time for this model (right panel). The times plotted on the left correspond to maximum shock extensions seen on 
the right. Here, the ante-sonic condition is computed as (c^/v'^g^)/ (0.187c). Instability should occur when this quantity exceeds approximately one (gray line in 
left panel). As can be seen in this figure, this occurs only during the final expansion of the shock (black line). 



finer resolution led to earlier explosions for the same driving 
neutrino luminosity. We have examined the dependence of 
both our ID and 2D results on resolution and show the results 
in Figure 12. While the development of small-scale instabili- 
ties and turbulence will depend heavily on the resolution used, 
for the present study the important question is whether or not 
changing the resolution significantly alter the time at which a 
model explodes for a given neutrino luminosity. To determine 
this, we ran several models using STOS with an increased 
minimum grid spacing at the maximum level of refinement of 
0.7 km, versus 0.5 km for our fiducial models. We find the 
results for explosion times as a function of neutrino luminosity 
are essentially identical in ID and not significantly different 
in 2D. As shown in Figure 12, the ID curves lie on top of one 
another and the 2D curves are very similar. For 2D, the small 
differences are that at lower driving luminosities (and, hence, 
later explosion times) the low-resolution models explode ear- 
lier, while at higher luminosities the lower resolution models 
explode later. Thus the high-resolution and low-resolution 
curves cross each other in Figure 12. 

It is difficult to make a direct comparison between our re- 
sults and those of Hanke et al. (201 1) because the nature of the 
computational grids employed is quite different. Hanke et al. 
(201 1) use 2D spherical-polar coordinates with non-equidistant 
spacing in radius. Here we use 2D cylindrical R — z coordi- 
nates with adaptive mesh refinement.^ Hanke et al. use a 
fiducial resolution of 3° in angle and 400 zones in radius. At a 
radius of 100 km, near the base of the gain region, this yields 
a linear grid spacing in the angular direction of 1.67 km, more 
than three times our fiducial resolution. For a spherical grid, 
this grid spacing will increase with radius while in our sim- 
ulations we allow refinement to the highest level anywhere 
in the computational domain, allowing us to achieve 0.5 km 
resolution for the entire post-shock region at all times. The 
highest resolution simulations of Hanke et al. use 0.5° angular 
resolution, resulting in 0.28 km linear resolution at 100 km and 
1 . 1 km linear resolution at 400 km. This resolution is likely 
comparable to our fiducial resolution, the linear grid spacing 
being smaller at the base of the gain region and larger at radii 
where shock oscillations reach the greatest extent. 

^ Nordhaus et al. (2010) also use 2D cylindrical coordinates with AMR. 
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Figure 12. Neutrino luminosity versus mass accretion rate for ID and 2D 
simulations with STOS at two different resolutions: the fiducial 0.5 km and 
0.7 km. For ID, the dash-dot curve represents using the first time the shock 
crosses 400 km as the explosion time for the L^^^^2 =2.4 model. Doing this 
makes the curves at different resolutions essentially identical. In both the ID 
and 2D cases, the resulting curves shown here are very similar, indicating our 
results are roughly converged in spatial resolution. 



We should not expect hydrodynamic simulations that involve 
instabilities and turbulence, such as these, to converge in spatial 
resolution, but the fact that we do not see significant differences 
in the explosion times at different resolutions indicates that our 
results are converged with respect to the result we seek. 

4. CONCLUSIONS 

We have conducted a parameter study exploring the depen- 
dence of the neutrino mechanism of core-collapse supemovae 
on the equation of state. We varied the driving neutrino lumi- 
nosity and EOS in several ID and 2D simulations of stellar 
core collapse and measured the resulting explosion times and 
mass accretion rates at the time of explosion. Table 2 and Fig- 
ure 3 summarize our results. We find that for every EOS, 2D 
explosion are obtained more easily than in ID, and that models 
using the Lattimer & Swesty EOS explodes more easily than 
models using the EOS of Shen et al., with LSI 80 resulting 
in slightly earlier explosions than LS220. Thus, our results 
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show a general trend that softer EOS lead to easier explosions. 
In order to test if the incompressibility parameter, K, is the 
primary determinant of the time delay to explosion for a given 
neutrino luminosity, we also ran a series of ID simulations 
using LS375. This model has a substantially higher incom- 
pressibility than STOS (375 MeV v. 281 MeV). We find that 
STOS, despite have a lower incompressibility, explodes later 
than LS375 for all neutrino luminosities. This highlights that 
the higher-order EOS parameters, such as symmetry energy 
coefficient, play an important role in determining the stiffness 
of an EOS. Our results of explosion delay times do follow a 
trend in resulting neutron star radii (see Fig. 1). EOS that 
produce larger typical neutron star radii explode later in our 
simulations. The dependence of the time-to-explosion on the 
EOS is likely due to the dependence of the acoustic energy flux 
from the PNS on the stiffness of the EOS. The softer EOS of 
Lattimer & Swesty result in more compact PNS that are more 
effective at translating advected perturbations into outgoing 
sound waves. This then results in more acoustic excitation of 
the standing accretion shock instability and an overall faster 
rate of shock expansion. This should be true no matter if the 
underlying SASI mechanism is acoustic (Blondin & Mezza- 
cappa 2006) or advective-acoustic (Foglizzo et al. 2007) since 
both mechanisms depend on the acoustic flux from the PNS. 
For high-L^^ models that do not depend strongly on the SASI, 
we attribute the difference between STOS and LS to the higher 
mass accretion rate just after bounce for STOS. 

In this study, we use a parameterized neutrino "light-bulb" 
with a fixed neutrino luminosity. There is no back-reaction 
from the EOS on the neutrino luminosity, allowing us to ex- 
amine the dependence of the resulting explosion times on the 
EOS alone. Previous studies have shown that the neutrino lumi- 
nosity can be significantly dependent on the EOS (Thompson 
et al. 2003; Sumiyoshi et al. 2005) and this is a very important 
consequence of the choice of EOS, but our study shows that 
even for a fixed neutrino luminosity, the likelihood of a given 
model to explode also depends on the EOS. 

The EOS for matter at temperatures and densities relevant 
to CCSNe is still an active area of research. As summarized 
in the Introduction and Table 1, available experiments and 
observations that constrain the EOS parameters favor LS220. 
LSI 80 is ruled out as it cannot produce a ~ 2Mq neutron star 
as required by the observations of PSR Jl 6 14-2230 (Demorest 
et al. 2010). STOS yields neutron stars of typical mass (around 
IAMq) with too large a radius (~ 15 km) to account for the 
results of Steiner et al. (2010) who estimate typical NS radii 
of 11-12 km. Several of the EOS parameters listed in Table 
1 for STOS also fall outside of the experimental and theoret- 
ical estimates. So, LS220 is the EOS that satisfies, or nearly 
satisfies, the largest number of constraints. Our results show, 
additionally, that LS220 is intermediate between STOS and 
LS 1 80 in favoring explosion for a given neutrino luminosity. 
Until we are more certain about the EOS at temperatures and 
densities relevant to CCSN simulations, LS220 seems a good 
choice. Some recent works have chosen to use LS220 on these 
same grounds (Ott et al. 2012; Dessart et al. 2012). 
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